home *** CD-ROM | disk | FTP | other *** search
/ CD Actual 3 / CD ACTUAL 3.iso / linux / incoming / libgr-2.000 / libgr-2 / libgr-2.0.3 / rle / hilbert.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-08-14  |  9.2 KB  |  411 lines
  1. /*
  2.  * This software is copyrighted as noted below.  It may be freely copied,
  3.  * modified, and redistributed, provided that the copyright notice is 
  4.  * preserved on all copies.
  5.  * 
  6.  * There is no warranty or other guarantee of fitness for this software,
  7.  * it is provided solely "as is".  Bug reports or fixes may be sent
  8.  * to the author, who may or may not act on them as he desires.
  9.  *
  10.  * You may not include this software in a program or other software product
  11.  * without supplying the source, or without informing the end-user that the 
  12.  * source is available for no extra charge.
  13.  *
  14.  * If you modify this software, you should include a notice giving the
  15.  * name of the person performing the modification, the date of modification,
  16.  * and the reason for such modification.
  17.  */
  18. /* 
  19.  * hilbert.c - Computes Hilbert curve coordinates from position and v.v.
  20.  * 
  21.  * Author:    Spencer W. Thomas
  22.  *         EECS Dept.
  23.  *         University of Michigan
  24.  * Date:    Thu Feb  7 1991
  25.  * Copyright (c) 1991, University of Michigan
  26.  */
  27. static char rcsid[] = "/usr/local/src/CVS/libgr/rle/hilbert.c,v 1.1.1.1 1995/08/15 02:46:22 neal Exp";
  28.  
  29. /*
  30.  * Lots of tables to simplify calculations.  Notation: p#i means bit i
  31.  * in byte p (high order bit first).
  32.  * p_to_s:    Output s is a byte from input p such that
  33.  *         s#i = p#i xor p#(i-1)
  34.  * s_to_p:    The inverse of the above.
  35.  * p_to_J:    Output J is "principle position" of input p.  The
  36.  *         principle position is the last bit s.t.
  37.  *         p#J != p#(n-1) (or n-1 if all bits are equal).
  38.  * bit:        bit[i] == (1 << (n - i))
  39.  * circshift:    circshift[b][i] is a right circular shift of b by i
  40.  *         bits in n bits.
  41.  * parity:    Parity[i] is 1 or 0 depending on the parity of i (1 is odd).
  42.  * bitof:    bitof[b][i] is b#i.
  43.  * nbits:    The value of n for which the above tables have been
  44.  *         calculated.
  45.  */
  46.  
  47. typedef unsigned int byte;
  48.  
  49. static int nbits = 0;
  50.  
  51. static byte
  52.     p_to_s[512],
  53.     s_to_p[512],
  54.     p_to_J[512],
  55.     bit[9],
  56.     circshift[512][9],
  57.     parity[512],
  58.     bitof[512][9];
  59.  
  60. /* Calculate the above tables when nbits changes. */
  61. static void
  62. calctables(n)
  63. int n;
  64. {
  65.     register int i, b;
  66.     int two_n = 1 << n;
  67.  
  68.     if ( nbits == n )
  69.     return;
  70.  
  71.     nbits = n;
  72.     /* bit array is easy. */
  73.     for ( b = 0; b < n; b++ )
  74.     bit[b] = 1 << (n - b - 1);
  75.  
  76.     /* Next, do bitof. */
  77.     for ( i = 0; i < two_n; i++ )
  78.     for ( b = 0; b < n; b++ )
  79.         bitof[i][b] = (i & bit[b]) ? 1 : 0;
  80.  
  81.     /* circshift is independent of the others. */
  82.     for ( i = 0; i < two_n; i++ )
  83.     for ( b = 0; b < n; b++ )
  84.         circshift[i][b] = (i >> (b)) |
  85.         ((i << (n - b)) & (two_n - 1));
  86.  
  87.     /* So is parity. */
  88.     parity[0] = 0;
  89.     for ( i = 1, b = 1; i < two_n; i++ )
  90.     {
  91.     /* Parity of i is opposite of the number you get when you
  92.      * knock the high order bit off of i.
  93.      */
  94.     if ( i == b * 2 )
  95.         b *= 2;
  96.     parity[i] = !parity[i - b];
  97.     }
  98.  
  99.     /* Now do p_to_s, s_to_p, and p_to_J. */
  100.     for ( i = 0; i < two_n; i++ )
  101.     {
  102.     int s;
  103.  
  104.     s = i & bit[0];
  105.     for ( b = 1; b < n; b++ )
  106.         if ( bitof[i][b] ^ bitof[i][b-1] )
  107.         s |= bit[b];
  108.     p_to_s[i] = s;
  109.     s_to_p[s] = i;
  110.  
  111.     p_to_J[i] = n - 1;
  112.     for ( b = 0; b < n; b++ )
  113.         if ( bitof[i][b] != bitof[i][n-1] )
  114.         p_to_J[i] = b;
  115.     }
  116. }
  117.  
  118. /*****************************************************************
  119.  * TAG( hilbert_i2c )
  120.  * 
  121.  * Convert an index into a Hilbert curve to a set of coordinates.
  122.  * Inputs:
  123.  *     n:    Number of coordinate axes.
  124.  *     m:    Number of bits per axis.
  125.  *     r:    The index, contains n*m bits (so n*m must be <= 32).
  126.  * Outputs:
  127.  *     a:    The list of n coordinates, each with m bits.
  128.  * Assumptions:
  129.  *     n*m < (sizeof r) * (bits_per_byte), n <= 8, m <= 9.
  130.  * Algorithm:
  131.  *     From A. R. Butz, "Alternative Algorithm for Hilbert's
  132.  *         Space-Filling Curve", IEEE Trans. Comp., April, 1971,
  133.  *         pp 424-426.
  134.  */
  135. void
  136. hilbert_i2c( n, m, r, a)
  137. int n;
  138. int m;
  139. long int r;
  140. int a[];
  141. {
  142.     byte rho[9], rh, J, sigma, tau,
  143.          sigmaT, tauT, tauT1 = 0, omega, omega1 = 0, alpha[9];
  144.     register int i, b;
  145.     int Jsum;
  146.  
  147.     /* Initialize bit twiddle tables. */
  148.     calctables(n);
  149.  
  150.     /* Distribute bits from r into rho. */
  151.     for ( i = m - 1; i >= 0; i-- )
  152.     {
  153.     rho[i] = r & ((1 << n) - 1);
  154.     r >>= n;
  155.     }
  156.  
  157.     /* Loop over bytes. */
  158.     Jsum = 0;
  159.     for ( i = 0; i < m; i++ )
  160.     {
  161.     rh = rho[i];
  162.     /* J[i] is principle position of rho[i]. */
  163.     J = p_to_J[rh];
  164.  
  165.     /* sigma[i] is derived from rho[i] by exclusive-oring adjacent bits. */
  166.     sigma = p_to_s[rh];
  167.  
  168.     /* tau[i] complements low bit of sigma[i], and bit at J[i] if
  169.      * necessary to make even parity.
  170.      */
  171.     tau = sigma ^ 1;
  172.     if ( parity[tau] )
  173.         tau ^= bit[J];
  174.     
  175.     /* sigmaT[i] is circular shift of sigma[i] by sum of J[0..i-1] */
  176.     /* tauT[i] is same circular shift of tau[i]. */
  177.     if ( Jsum > 0 )
  178.     {
  179.         sigmaT = circshift[sigma][Jsum];
  180.         tauT = circshift[tau][Jsum];
  181.     }
  182.     else
  183.     {
  184.         sigmaT = sigma;
  185.         tauT = tau;
  186.     }
  187.  
  188.     Jsum += J;
  189.     if ( Jsum >= n )
  190.         Jsum -= n;
  191.  
  192.     /* omega[i] is xor of omega[i-1] and tauT[i-1]. */
  193.     if ( i == 0 )
  194.         omega = 0;
  195.     else
  196.         omega = omega1 ^ tauT1;
  197.     omega1 = omega;
  198.     tauT1 = tauT;
  199.  
  200.     /* alpha[i] is xor of omega[i] and sigmaT[i] */
  201.     alpha[i] = omega ^ sigmaT;
  202.     }
  203.  
  204.     /* Build coordinates by taking bits from alphas. */
  205.     for ( b = 0; b < n; b++ )
  206.     {
  207.     register int ab, bt;
  208.     ab = 0;
  209.     bt = bit[b];
  210.     /* Unroll the loop that stuffs bits into ab.
  211.      * The result is shifted left by 9-m bits.
  212.      */
  213.     switch( m )
  214.     {
  215.     case 9:    if ( alpha[8] & bt) ab |= 0x01;
  216.     case 8:    if ( alpha[7] & bt) ab |= 0x02;
  217.     case 7:    if ( alpha[6] & bt) ab |= 0x04;
  218.     case 6:    if ( alpha[5] & bt) ab |= 0x08;
  219.     case 5:    if ( alpha[4] & bt) ab |= 0x10;
  220.     case 4:    if ( alpha[3] & bt) ab |= 0x20;
  221.     case 3:    if ( alpha[2] & bt) ab |= 0x40;
  222.     case 2:    if ( alpha[1] & bt) ab |= 0x80;
  223.     case 1:    if ( alpha[0] & bt) ab |= 0x100;
  224.     }
  225.     a[b] = ab >> (9 - m);
  226.     }
  227. }
  228.  
  229. /*****************************************************************
  230.  * TAG( hilbert_c2i )
  231.  * 
  232.  * Convert coordinates of a point on a Hilbert curve to its index.
  233.  * Inputs:
  234.  *     n:    Number of coordinates.
  235.  *     m:    Number of bits/coordinate.
  236.  *     a:    Array of n m-bit coordinates.
  237.  * Outputs:
  238.  *     r:    Output index value.  n*m bits.
  239.  * Assumptions:
  240.  *     n*m <= 32, n <= 8, m <= 9.
  241.  * Algorithm:
  242.  *     Invert the above.
  243.  */
  244. void
  245. hilbert_c2i( n, m, a, r)
  246. int n;
  247. int m;
  248. int a[];
  249. long int *r;
  250. {
  251.     byte rho[9], J, sigma, tau,
  252.          sigmaT, tauT, tauT1 = 0, omega, omega1 = 0, alpha[9];
  253.     register int i, b;
  254.     int Jsum;
  255.     long int rl;
  256.  
  257.     calctables(n);
  258.  
  259.     /* Unpack the coordinates into alpha[i]. */
  260.     /* First, zero out the alphas. */
  261.     switch ( m ) {
  262.     case 9: alpha[8] = 0;
  263.     case 8: alpha[7] = 0;
  264.     case 7: alpha[6] = 0;
  265.     case 6: alpha[5] = 0;
  266.     case 5: alpha[4] = 0;
  267.     case 4: alpha[3] = 0;
  268.     case 3: alpha[2] = 0;
  269.     case 2: alpha[1] = 0;
  270.     case 1: alpha[0] = 0;
  271.     }
  272.  
  273.     /* The loop that unpacks bits of a[b] into alpha[i] has been unrolled.
  274.      * The high-order bit of a[b] has to go into alpha[0], so we pre-shift
  275.      * a[b] so that its high-order bit is always in the 0x100 position.
  276.      */
  277.     for ( b = 0; b < n; b++ )
  278.     {
  279.     register int bt = bit[b], t = a[b] << (9 - m);
  280.  
  281.     switch (m)
  282.     {
  283.     case 9: if ( t & 0x01 ) alpha[8] |= bt;
  284.     case 8: if ( t & 0x02 ) alpha[7] |= bt;
  285.     case 7: if ( t & 0x04 ) alpha[6] |= bt;
  286.     case 6: if ( t & 0x08 ) alpha[5] |= bt;
  287.     case 5: if ( t & 0x10 ) alpha[4] |= bt;
  288.     case 4: if ( t & 0x20 ) alpha[3] |= bt;
  289.     case 3: if ( t & 0x40 ) alpha[2] |= bt;
  290.     case 2: if ( t & 0x80 ) alpha[1] |= bt;
  291.     case 1: if ( t & 0x100 ) alpha[0] |= bt;
  292.     }
  293.     }
  294.  
  295.     Jsum = 0;
  296.     for ( i = 0; i < m; i++ )
  297.     {
  298.     /* Compute omega[i] = omega[i-1] xor tauT[i-1]. */
  299.     if ( i == 0 )
  300.         omega = 0;
  301.     else
  302.         omega = omega1 ^ tauT1;
  303.  
  304.     sigmaT = alpha[i] ^ omega;
  305.     /* sigma[i] is the left circular shift of sigmaT[i]. */
  306.     if ( Jsum != 0 )
  307.         sigma = circshift[sigmaT][n - Jsum];
  308.     else
  309.         sigma = sigmaT;
  310.  
  311.     rho[i] = s_to_p[sigma];
  312.  
  313.     /* Now we can get the principle position. */
  314.     J = p_to_J[rho[i]];
  315.  
  316.     /* And compute tau[i] and tauT[i]. */
  317.     /* tau[i] complements low bit of sigma[i], and bit at J[i] if
  318.      * necessary to make even parity.
  319.      */
  320.     tau = sigma ^ 1;
  321.     if ( parity[tau] )
  322.         tau ^= bit[J];
  323.  
  324.     /* tauT[i] is right circular shift of tau[i]. */
  325.     if ( Jsum != 0 )
  326.         tauT = circshift[tau][Jsum];
  327.     else
  328.         tauT = tau;
  329.     Jsum += J;
  330.     if ( Jsum >= n )
  331.         Jsum -= n;
  332.  
  333.     /* Carry forth the "i-1" values. */
  334.     tauT1 = tauT;
  335.     omega1 = omega;
  336.     }
  337.  
  338.     /* Pack rho values into r. */
  339.     rl = 0;
  340.     for ( i = 0; i < m; i++ )
  341.     rl = (rl << n) | rho[i];
  342.     *r = rl;
  343. }
  344.  
  345.  
  346. #ifdef test
  347. #include <stdio.h>
  348.  
  349. main()
  350. {
  351.     int a[9], n, m, i;
  352.     long int r, r1;
  353.  
  354.     for (;;)
  355.     {
  356.     printf( "Enter n, m: " );
  357.     scanf( "%d", &n );
  358.     if ( n == 0 )
  359.         break;
  360.     scanf( "%d", &m );
  361.     while ( (i = getchar()) != '\n' && i != EOF )
  362.         ;
  363.     if ( i == EOF )
  364.         break;
  365.     for ( r = 0; r < 1 << (n*m); r++ )
  366.     {
  367.         hilbert_i2c( n, m, r, a );
  368.         if ( n*m <= 6 )
  369.         {
  370.         printf( "a = " );
  371.         for ( i = 0; i < n; i++ )
  372.             printf( "0x%0*x ", (m+3)/4, a[i] );
  373.         }
  374.         hilbert_c2i( n, m, a, &r1 );
  375.         if ( n*m <= 6 )
  376.         printf( "r = 0x%0*x\n", (n*m+3)/4, r1 );
  377.         else if ( r != r1 )
  378.         printf( "r = 0x%0*x; r1 = 0x%0*x\n", (n*m+3)/4, r,
  379.             (n*m+3)/4, r1 );
  380.     }
  381.     }
  382. }
  383.  
  384. p1t( tbl, n )
  385. byte tbl[];
  386. int n;
  387. {
  388.     int i;
  389.  
  390.     for ( i = 0; i < n; i++ )
  391.     printf( "%02x ", tbl[i] );
  392.     putchar( '\n' );
  393. }
  394.  
  395. p2t( tbl, n, m )
  396. byte tbl[][9];
  397. int n;
  398. {
  399.     int i;
  400.  
  401.     for ( i = 0; i < n; i++ )
  402.     {
  403.     printf( "%3d: ", i );
  404.     p1t( tbl[i], m );
  405.     }
  406. }
  407.  
  408.  
  409. #endif
  410.  
  411.